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Abstract 

We  consider  the  use  of  a  multigrid  method  with  central  differencing  to  solve  the 
Navier-Stokes  equations  for  hypersonic  flows.  The  time-dependent  form  of  the  equations 
is  integrated  with  an  explicit  Runge-Kutta  scheme  accelerated  by  local  time  stepping  and 
implicit  residual  smoothing.  Variable  coefficients  are  developed  for  the  implicit  process 
that  remove  the  diffusion  limit  on  the  time  step,  producing  significant  improvement  in 
convergence.  A  numerical  dissipation  formulation  that  provides  good  shock-capturing 
capability  for  hypersonic  flows  is  presented.  This  formulation  is  shown  to  be  a  crucial 
aspect  of  the  multigrid  method.  Solutions  are  given  for  two-dimensional  viscous  flow 
over  a  NACA  0012  airfoil  and  three-dimensional  viscous  flow  over  a  blunt  biconic. 

’This  research  was  partially  supported  by  the  National  Aeronautics  and  Space  Administration  under 
NASA  Contract  No.  NASl-18605  while  the  first  author  was  in  residence  at  the  Institute  for  Computer 
Applications  in  Science  and  Engineering  (ICASE)  NASA  Langley  Research  Center,  Hampton,  VA  23665. 
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Introduction 

At  the  present  time  there  is  a  strong  interest  in  high-speed  flight  vehicles.  Some 
examples  of  these  vehicles  are  the  high-speed  civil  transport  (HSCT)  and  the  hypersonic 
flight  configurations  being  considered  for  me  National  Aero-Space  Plane  (NASP).  In  the 
Ccise  of  the  NASP,  one  encounters  complex  high  Mach  number  phenomena  and  interac¬ 
tions,  which  can  involve  strong  shock  and  expansion  waves,  in  not  only  the  flow  over 
the  vehicle,  but  also  in  the  flow  through  the  engines,  where  chemical  reactions  occur. 
An  effective  design  method  for  such  vehicles  will  obviously  require  both  detailed  experi¬ 
mental  data  as  well  as  flexible,  efficient,  and  accurate  computational  techniques.  Robust 
prediction  methods  (i.e.,  those  that  con  be  applied  on  a  routine  basis)  are  not  currently 
available  for  hypersonic  flows.  However,  due  to  the  significant  progress  during  the  last 
decade  in  the  development  of  effective  algorithms  for  subsonic  and  transonic  flows,  there 
are  a  number  of  oppoitunities  for  constructing  improved  schemes  for  high-speed  flows. 

One  powerful  approach  for  the  numerical  solution  of  partial  differential  equations, 
which  has  been  successfully  applied  to  fluid  flow  problems,  is  multigrid.  Multigrid  meth¬ 
ods  were  first  developed  for  elliptic  equations.  These  were  later  extended  to  hyperbolic 
equations  such  as  the  time-dependent  fluid  dynamic  equations  for  subsonic  and  tran¬ 
sonic  flow  [1,  2,  3,  4,  5].  Even  for  transonic  crises,  the  steady  state  can  retain  many  of 
the  properties  of  an  elliptic  equation  when  the  region  of  supersonic  flow  is  limited.  We 
shall  show  that  with  proper  care  the  multigrid  method  still  works  for  hypersonic  flow. 
Gustafsson  and  Lotstedt  [6]  have  pointed  out  that  hyperbolic  multigrid  works  by  two 
different  processes.  For  the  long  waves,  the  advection  process  is  most  important  and 
multigrid  achieves  its  efficiency  by  allowing  the  use  of  larger  time  steps  on  coarser  grids. 
Hence,  it  is  important  that  the  smoother  use  large  time  steps.  However,  for  the  shorter 
waves,  dissipation  is  more  important  and  the  efficiency  of  multigrid  is  based  on  principles 
similar  to  that  for  elliptic  equations. 

Several  investigators  have  applied  multigrid  to  high-speed  flows  with  varying  degrees 
of  success.  For  example,  in  [7]  the  Euler  equations  are  solved  with  and  without  chemistry 
using  the  Roe  scheme  for  spatial  discretization  and  the  ADI  scheme  for  time  marching. 
A  factor  of  four  decrease  in  computational  time  was  obtained  with  the  multigrid  method 
for  a  simple  one-dimensional  nozzle  flow  (exit  Mach  number  of  about  3.7).  A  calculation 
for  a  Mach  22  wedge  flow  showed  the  basic  scheme  to  be  noticeably  faster  than  the 
multigrid  scheme.  One  notable  conclusion  of  this  work  was  that  the  performance  of  the 
multigrid  is  driven  by  the  fluid  dynamics  and  not  the  chemistry  (at  least  foi  the  case  of 
simple  reaction  models).  The  high  level  of  performance  and  widespread  application  of 
multigrid  algorithms  with  central  diffwoncing  and  an  explicit  multistage  time-stepping 
scheme  have  provided  strong  encouragement  to  make  them  work  for  hypersonic  problems. 
An  initial  effort  [8]  to  apply  this  type  of  algorithm  resulted  in  numerical  difficulties  that 
prevented  the  calculation  of  two-dimensional  flows  (i.e.,  blunt  body  and  wedge  type)  with 
a  Mach  number  higher  than  about  7.  In  order  to  compute  such  flows,  a  low  Courant- 
Friedrichs-Lewy  (CFL)  number  was  required.  Thus  four  and  five  stage  schemes  were  not 
practical,  since  there  is  substantial  deterioration  in  the  high  frequency  damping  of  the 
scheme  due  to  the  large  reduction  in  the  CFL  number.  The  CFL  restriction  reduced  the 
potential  of  the  scheme  as  a  viscous  flow  solver.  More  recently  an  algorithm  utilizing  a 
semicoarsening  technique,  a  symmetric  TVD  formulation,  and  a  three  stage  Runge-Kutta 
scheme  [9]  was  proposed  and  used  to  compute  high  Reynolds  number  (laminar)  Mach 
10  flow  over  an  airfoil  at  10  degrees  angle  of  attack.  A  good  resolution  of  the  bow  shock 


1 


wave  and  a  reasonable  convergence  rate  were  obtained.  The  method  of  semicoarsening 
considered  required  a  much  more  complicated  cycle  strategy  than  that  employed  with 
standard  multigrid  methods.  In  addition,  it  appears  to  be  somewhat  cumbersome  to 
implement  in  three  dimensions. 

It  is  our  contention  that  standard  multigrid  techniques  can  be  used  in  conjunction 
with  central  differencing  for  hypersonic  flows.  To  do  this  one  must  ensure  that  the 
basic  algorithm  exhibits  good  damping  of  high  frequencies,  on  both  the  fine  and  coarse 
meshes,  in  the  neighborhood  of  strong  discontinuities.  It  becomes  more  difficult  to 
eliminate  these  high  frequency  oscillations  as  the  Mach  number  increases,  since  the  jumps 
across  shocks  become  larger.  Thus  a  considerable  part  of  the  following  discussion  will 
concentrate  on  the  smoothing  algorithm.  The  fundamental  features  of  the  multigrid 
process  (i.e..  Full  Approximation  Storage  scheme,  grid  transfer  operators,  fixed  cycle 
strategy)  are  fairly  standard.  Other  aspects,  such  as  type  of  coarse  grid  correction 
scheme  and  procedure  for  smoothing  of  coarse  grid  corrections,  found  crucial  in  the 
present  work  will  be  emphasized.  In  this  paper  we  consider  a  Runge-Kutta  scheme  [10] 
<.s  the  smoother  foi  the  multigrid  method.  Central  differences  for  spatial  approximations 
are  augnicnted  by  an  artificial  viscosity  bcised  on  TVD  principles  [11].  Several  changes 
are  made  to  the  numerical  algorithm  so  that  a  converged  solution  can  be  obtained  for 
high-speed  flows.  Initially,  we  describe  the  bcisic  Runge-Kutta  method  for  the  central- 
difference  scheme  with  the  numerical  viscosity.  We  discuss  an  implicit  procedure  for 
accelerating  the  convergence  of  the  basic  scheme.  This  procedure  removes  the  diffusion 
limit  on  the  time  step,  which  can  be  quite  severe  at  hypersonic  speeds.  We  then  apply 
the  multigrid  method  to  solve  both  two-  and  three-dimensional  viscous  flows. 

Basic  Scheme 

The  basic  elements  vf  the  scalar  dissipation  model  considered  in  this  paper  wer 
first  introduced  by  Jameson,  Schmidt,  and  Turkel  [10]  in  conjunction  with  Runge-Kutta 
explicit  schemes.  The  spatial  discretization  is  based  on  central  differences  with  an  addi¬ 
tional  artificial  viscosity.  This  algorithm  has  been  used  by  many  investigators  to  solve  the 
Euler  equations  numerically  for  a  wide  range  of  fluid  dynamic  applications.  The  same 
type  of  spatial  discretization  hcis  been  applied  to  alternating  direction  implicit  (ADI) 
schemes  [12]  and  LU  factored  implicit  schemes  [13].  In  this  section  the  basic  scheme 
is  briefly  reviewed.  Since  the  elements  of  the  two-dimensional  and  three-dimensional 
schemes  are  essentially  the  same,  we  will  take,  for  the  purpose  of  simplicity,  a  two- 
dimensional  point  of  view  in  describing  the  basic  scheme. 

Consider  the  Euler  equations  in  the  form 

+  /x  +  <7y  =  0>  (1) 

where  the  four-component  vector  of  conserved  variables 

l'K  =  [.p  pu  pv  pEf,  (2) 

and  f,g  are  the  corresponding  flux  vectors.  The  quantity  p  is  the  density,  u  and  v 
are  the  Cartesian  velocity  components,  and  E  is  the  sf  ecific  total  internal  energy.  The 
independent  variables  are  time  t  and  Cartesian  coordinates  (x,  y).  If  (1)  is  transformed 
to  arbitrary  curvilinear  coordinates  ^  =  {(.x,?/)  and  t}  =  i]{x,y),.then  we  obtain 

{J-^W),  +  F^  +  G„  =  0  (3) 
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where  J~^  is  the  inverse  transformation  Jacobian,  and 

F  =  fVn-  G  =  gx^-  fy^. 

In  a  cell-centered,  finite-volume  method,  (1)  is  integrated  over  an  elemental  volume  in 
the  discretized  computational  domain,  and  J~^  is  identified  as  the  volume  of  the  cell. 
Equation  (3)  can  also  be  written  as 

J-^Wt  -i-  AW^  +  BWr,  =  0 


where  A  and  B  are  the  flux  Jacobian  matrices  defined  by  A  =  dFjdW  and  B  =  dGIdW. 

To  advance  the  scheme  in  time  we  use  a  multistage  scheme.  A  typical  step  of  a 
Runge-Kutta  approximation  to  (3)  is 

-  AD\  ,  (4) 

where  and  Dj,  are  spatial  differencing  operators,  and  AD  represents  the  artificial 
dissipation  terms.  The  derivatives  of  the  fluxes  are  approximated  by  central  differences. 
Hence,  we  need  to  evaluate  F  at  {i  +  j).  For  the  physical  diffusive  fluxes  the  approxima¬ 
tions  follow  directly.  Since  the  dependent  variables  are  given  at  {i,j)  we  need  to  average 
them  to  calculate  the  convective  flux  at  the  cell  face.  This  averaging  can  be  done  in 
several  ways.  In  the  original  scheme  the  conservative  variables  were  averaged  to  find  the 
variables  at  the  cell  face  and  then  the  flux  was  evaluated.  Moreover,  in  two  dimensions 
five  quantities  were  averaged,  the  four  conservative  variables  and  pE +p  .  This  was  done 
because  enthalpy  damping  was  used.  Alternatively  one  can  average  density,  momentum, 
and  total  enthalpy.  This  might  be  more  accurate  for  the  Euler  equations  since  the  total 
enthalpy  is  constant  in  the  steady  state.  An  alternative  to  this  procedure  is  to  average 
the  fluxes  rather  than  the  dependent  variables,  where  the  fluxes  are  defined  as  for  (3)  and 
the  transformation  derivatives  are  evaluated  at  the  cell  interface.  Averaging  the  fluxes 
allows  for  a  more  accurate  representation  of  the  Rankine-Hugoniot  jump  conditions  in 
one  dimension.  In  this  paper  we  average  the  fluxes  but  have  not  seen  much  difference 
between  the  various  approaches. 

The  dissipation  terms  are  a  blending  of  second  and  fourth  differences.  That  is, 


where 


AD  =  {dI  +  Dl~  Dl  -  Di)  W, 

(5) 

(6) 

Z)|1V  =  Vf  AfVfA;]  W^j, 

(7) 

and  are  the  standard  forward  and  backward  difference  operators,  respectively, 

associated  with  the  ^  direction.  The  variable  scaling  factor  A  is  chosen  as 


where  A^  is  proportional  to  the  spectral  radius  of  the  matrix  A.  The  coefficients  and 
are  adapted  to  the  flow  and  are  defined  as  follows: 


J2)  ^ 


max(;/._,,j, 


(9) 
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= 


Pi+i,i  -  ^Pi,j  +  Pi-1, j 


(-i) 

e.,  1  .  =  max 

»+2.J 


Pi+i,i  +  ‘^Pi,]  +  Pi-1, i  1  ’ 


(10) 

(11) 


where  p  is  the  pressure,  and  the  quantities  and  are  constants  to  be  specified. 
The  operators  for  the  tj  direction  are  defined  in  a  similar  manner. 

The  second-difference  dissipation  term  is  nonlinear.  Its  purpose  is  to  introduce  an 
entropy-like  condition  and  to  suppress  oscillations  in  the  neighborhood  of  shocks.  This 
term  is  small  in  the  smooth  portion  of  the  flow  field.  The  fourth-difference  dissipation 
term  is  basically  linear  and  is  included  to  damp  high-frequency  modes  and  allow  the 
scheme  to  approach  a  steady  state.  Only  this  term  affects  the  linear  stability  of  the 
scheme.  Near  shocks  it  is  reduced  to  zero.  For  high-speed  flows,  the  switch  (30)  is  not 
very  good  and  does  not  allow  the  multigrid  to  converge.  Instead  we  consider  a  TVD 
variation  of  the  switch  [11]  given  by 


,  _  \Pi+i,j  ^Pi,j  ~b  _  (2)  _  1  /9 

\Pi^i,i  -  Pi, i\  +  \Pi,i- Pi-1, i\  +  ^^ 


(12) 


With  this  change  and  the  factor  1/2  in  front  of  the  second-difference  dissipation  term, 
the  scalar  equation  becomes  first-order  upwind  near  shocks.  In  the  case  of  the  original 
I/,  we  find  that  u  ~  0.05  near  shock  waves  in  transonic  flows.  The  parameter  e  must  be 
chosen  carefully  to  prevent  the  switch  from  being  activated  by  noise.  In  fact,  we  found 
it  useful  to  take  an  average  of  the  two  versions  for  v.  Hence,  we  use 


_  \Pi+i,3-‘^Pi,i  +  Pi-i,i\ 

(1  -  w)  {PTVD)i,j  +  ’ 

where 


(13) 


{PTVD)i,i  —  lpi+l,i  Pi,i\  "b  \pi,j 
Pi,j  =  Pm,)  +  2p;.j  -b  Pi-I,j, 

and  w  is  taken  to  be  1/2.  We  now  no  longer  have  a  free  parameter  for  the  second- 
difference  dissipation. 

Several  other  changes  were  made  to  the  scheme  in  addition  to  the  change  to  a  TVD 
switch.  In  the  original  algorithm,  the  artificial  viscosity  for  the  energy  equation  was 
based  on  the  total  enthalpy  rather  than  the  total  internal  energy.  For  high-speed  flows 
we  base  the  artificial  viscosity  on  the  total  internal  energy  so  that  in  each  equation  the 
basic  dependent  variable  is  also  used  in  the  artificial  viscosity.  This  is  more  in  line  with 
upwind  schemes.  This  has  previously  been  used  in  central-difference  schemes  [14].  The 
algorithm  no  longer  preserves  a  constant  total  enthalpy  in  the  steady  state  (eis  the  Euler 
equations  do),  but  enthalpy  damping  is  not  useful  for  supersonic  flows.  In  most  cases  the 
difference  between  the  two  approaches  is  small  with  each  approach  having  its  advantages. 
The  original  version  seems  to  give  slightly  sharper  shocks,  while  the  other  one  appears 
to  make  the  scheme  more  robust. 

The  form  of  the  dissipation  model  of  the  basic  or  driving  scheme  is  usually  modified 
for  the  coarse  grids  in  the  multigrid  process.  A  constant  coefficient,  second-difference  dis¬ 
sipation  is  not  only  less  expensive  computationally  but  also  generally  provides  adequate 
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smoothing  properti^.  For  high-speed  flows  we  find  it  necessary  to  append  a  rionjiriear 
dissipation  to  the  usual  one  that  depends  on  the  modified  switching  function  of  (12).  We 
also  need  to  increeise  the  constant  coefficient  on  a  coarse  mesh  from  the  standard  value 
of  1/16  to  a  value  of  1/4. 

In  order  for  the  scheme  to  be  stable  it  is  necessary  to  restrict  the  time  step.  Both  a 
convection  limit  as  well  as  a  diffusion  limit  must  be  taken  into  account  in  general.  As 
shown  in  [15]  the  actual  time  step  (Atact),  based  on  a  sufficient  condition  for  stability,  is 
determined  as  follows: 

Atact  <  +  A,  -b  dA„,jcoua]  j  (14) 

where  N  is  taken  to  be  the  allowable  CFL  number,  and  the  constant  d  is  4.  The  spectral 
radii  A^  and  A,,  of  the  Jacobian  matrices  A  and  B,  where  the  tilde  denotes  that  the 
matrix  is  multiplied  by  the  transformation  Jacobian  J.  are  given  by 

A4  =  +  C)/5T^, 


A^  =  Hy  +  UT)^\  +  c^r]l+T)l. 
For  the  thin-layer  Navier-Stokes  equations. 


7/t 

Prp 


ivl  +  vl). 


where  7  is  the  specific  heat  ratio,  and  Pr  is  the  Prandtl  number. 

We  now  consider  how  the  contributions  to  (14)  behave  near  the  body  as  the  inflow 
Mach  number  is  increased.  We  assume  that  the  variables  are  normalized  so  that  the 
density  and  pressure  are  1  at  inflow.  In  these  nondimensional  units  the  viscous  te-ms 
are  multiplied  by  In  the  boundary  layer  the  velocity  goes  to  zero  and  so  A, 

basically  behaves  as  the  sound  speed  c.  If  /i  =  T  =  p/p,  then 


VISCOUS 

TT 


p^RccoPr 


We  next  consider  the  ratios  pi^/pi  and  pt^/pi  for  adiabatic  flow  (16).  As  Moo  becomes 
larger  p  behaves  like  QM^jo  while  p  approaches  6.  So  as  Moo  becomes  larger 


viscous 

K 


Thus,  the  diffusion  limit  on  the  time  step  is  quite  significant  for  high  Mach  numbers. 

For  all  flow  calculations  in  this  paper  a  five  stage  Runge-Kutta  scheme  with  a  weighted 
evaluation,  as  detailed  in  [17],  of  the  dissipation  terms  on  the  first,  third,  and  fifth  stages 
is  used.  The  time  step  is  reduced  near  shocks  by  including  a  term  that  depends  on  Q 

The  reduction  is  constructed  so  that  there  is  a  CFL  number  of  1  when  1/  =  1.  It  serves  to  ^ 

reduce  the  magnitude  of  the  change  in  the  solution  near  the  shock  wave,  which  exhibits  ' 

strong  nonlinear  behavior. 
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Implicit  B.esidual  Smoothing 

A  mathematical  step  can  be  included  in  the  basic  time-stepping  scheme  to  extend 
the  local  stability.  This  procedure  is  called  implicit  smoothing  (or  averaging)  of  the 
residuals.  Previously,  variable  coefficients  for  this  implicit  process  have  been  introduced 
by  Martinelli  [2]  and  Swanson  and  Turkel  [18].  In  addition,  an  explicit  multidimensional 
derivation  of  such  coefficients  is  considered  in  [19].  These  variable  coefficients  have 
proven  to  be  quite  reliable  in  extending  the  stability  limit  of  the  Runge-Kutta  scheme  by 
a  factor  of  two.  However,  their  development  is  based  on  hyperbolic  considerations  only, 
’d  they  do  not  eliminate  the  need  for  a  diffusion  factor  in  the  time  step.  As  indicated 
in  the  preceding  section  and  we  emphasize  here,  this  factor  can  even  dominate  the  time 
step  determination  in  the  C2ise  of  hypersonic  flow.  In  this  section  we  will  briefly  discuss 
the  smoothing  coefficients  of  [18]  and  present  their  straightforward  extension  to  three 
dimensions.  Then  we  will  develop  smoothing  coefficients  that  allow  the  removal  of  the 
diffusion  limit,  resulting  in  a  At  which  depends  only  on  the  spectral  radii  of  the  flux 
.Jacobian  matrices. 

For  multidimensional  problems,  the  residual  smoothing  can  be  applied  in  the  form 


n(/-/5,V,A,) 


'Hi  '^«ii 


l  =  u,c 


where  the  residual  7^;^^  is  defined  by 


(15) 


KS?  =  +  -COM'S’  -  (16) 

and  computed  in  the  Runge-Kutta  stage  m.  The  quantity  is  the  total  artificial 

dissipation  at  stage  m,  and  is  the  final  residual  at  stage  m  after  the  sequence  of 
smoothings  involving  each  coordinate  direction.  The  difference  operators  Cc  and  Co 
are  associated  with  the  convection  and  physical  diffusion  terms.  To  derive  the  maximum 
stability  extension  for  the  hyperbolic  problem,  the  implicit  procedure  is  applied  after 
each  stage  of  the  Runge-Kutta  scheme.  The  coefficients  are  functions  of  the  spectral 
radii  A^,  A,,,  and  A,^.  In  two  dimensions  they  are  written  as  follows: 


Pi  = 

P.  = 


max 


max 


jy_  1 
N‘  1  -{-^r,,£ 

JL  1 

N‘  l+i/ir-l 


)’ 


(17) 

(18) 


where  the  ratio  7*,,^  =  Ajj/A^,  and  the  quantity  N/N’  is  the  ratio  of  the  CFL  number 
of  the  smoothed  scheme  to  that  of  the  basic  explicit  scheme  (usually  having  a  value  of 
2).  In  hypersonic  flow  applications  we  found  it  necessary  for  N“  to  be  3.25,  rather  than 
the  value  of  3.75  used  for  transonic  computations.  From  a  linear  stability  analysis,  the 
scheme  with  these  coefficients  is  stable  for  all  mesh  cell  aspect  ratios  when  the  parameter 
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‘iff  J§  0fl2S  ^dliV/iV*  is  sufficiently  l^ge.  The  practice  liimtaUon  on  tHeiGdur^fenu^er 
is  due  to  the  requirement  for  effective  high,  frequency /damping.  For  large  iV/iV*,  the 
high  frequency  damping  of  the  scheme  vanishes.  The  variable  coefficients?are  functions 
of  the  local  mesh  cell  aspect  ratio,  and  thus  the  smoothing  process  is  not  activated 
in  a  coordinate  direction  where  it  is  not  needed.  This  is  important  for  best  possible 
convergence. 

The  formulas  of  (17)  and  (IS)  can  be  easily  extended  to  three  dimensions.  Moreover, 
we  define 


/3,  =  max|i  (^*<)  -1  .c}.  >  = 


(19) 


where 


= 


1 


l  +  ^3-z?(^’,{  +  r«)’ 
1 


l+^^’3-D(^•„^+r(,)’ 


= 


1 


+rr  ) 


Once  again  r  represents  a  ratio  of  characteristic  speeds.  A  typical  value  for  is 
0.0625. 

We  now  snow  how  one  can  utilize  these  coefficients  as  the  basis  for  new  coefficients 
that  will  remove  the  diffusion  limit  on  the  time  step.  Consider  the  scalar  diffusion 
equation 

dw  d^w 


dt  ^  dy^ 


(20) 


II  we  approximate  the  spatial  derivative  of  (20)  with  a  central  difference  and  take  a 
Fourier  transform,  we  obtain 

Ai^  =  Zw\  (21) 

where  the  caret  indicates  a  transformed  quantity,  and  the  Fourier  symbol 


Z  =  — 2Arf(l  —  co$0) 

with  Xd  =  AtfiflXy^.  If  implicit  residual  smoothing  is  applied,  then 

_  -2Aj(l  -  cosO) 

1  +  2pdil  —  cosO) 

A  sufficient  condition  for  stability  is  given  by 

\Z\<Nd  for  all  0, 


(22) 

(23) 

(24) 


wh^re  A^J  is  the  diffusion  number  of  the  unsmoothed  scheme.  It  then  follows  that  the 
smoothing  coefficient  Pd  is  given  by 


f>.>\ 


n;  , 


Atd 

Ail 


-1 


(25) 
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where  Nj  =  4Aii.  This  implicit  smoothing  coefficient  has  a  linear  dependence  on  the  dif¬ 
fusion  number  ratio.  In  the  case  of  the  pure  hyperbolic  problem,  one  obtains  a  quadratic 
dependence  with  the  CFL  number  ratio.  Thus  we  must  determine  a  w'ay  to  combine 
tliese  results  to  obtain  a  coefficient  that  is  valid  for  an  equation  with  both  convection 
and  diffusion  effects. 

Suppose  w'e  consider  the  thin-layer  form  of  the  tw'o-dimensional  Navier-Stokes  equa¬ 
tions,  then  one  could  simply  use  the  smoothing  coefficient  of  (17)  in  the  streamwise-like 
(^)  direction.  One  would  anticipate  that  a  possible  formulation  for  the  normal  tj  direc¬ 
tion  would  depend  on  a  diffusion  type  ,3  near  the  surface  and  a  convection  type  when 
the  viscous  effects  are  no  longer  important.  Notice  the  dependency  in  (25)  on  the  ratio 
of  the  actual  At  to  the  At  of  the  basic  explicit  scheme.  So  we  define 


Aigg 


(26) 


Since  we  want  to  remove  the  diffusion  limit,  the  actual  At  must  depend  only  on  and 
A,.  Thus  we  rewrite  (26)  as 


mn  = 


1 

4 


^  (■^ti)n 

.^'a^  +  a. 


which,  we  can  approximate  simply  as 


(27) 


It  does  not  appear  to  matter  whether  one  uses  the  form  of  (27)  or  that  of  (28),  provided 
the  constant  is  defined  properly.  Here  we  have  elected  to  use  the  simpler  form  of  (28), 
which  has  also  been  considered  by  Radespicl  and  Kroll  (20).  From  numerical  experiments 
we  found  that  a  satisfactory  x’alue  for  Ci  is  10. 

The  variable  coefficient  of  (28)  cannot  generally  be  used  by  itself.  For  example,  in 
an  airfoil  flow,  {^d)r,  goes  to  zero  too  fast  at  the  leading  edge,  resulting  in  a  zero  v'alue 
in  the  inviscid  region.  We  have  overcome  this  difficulty  by  calculating  0,,  as 


=  max  {{0d)„ ,  {Pc)r,) ,  (29) 

w’here  (/3c)tj  is  defined  by  (18).  In  the  Results  section  of  this  paper  we  consider  Mach  10 
and  Mdch  20  turbulent  flow  over  an  airfoil.  The  removal  of  the  diffusion  limit  allowed 
the  residual  in  each  case  to  be  reduced  at  least  an  additional  order  of  magnitude  in  200 
multigrid  cycles. 


Multigrid  Method 

.4s  indicated  earlier  the  salient  features  of  the  multigrid  method  considered  here  are 
fairly  standard.  We  apply  the  Full  Approximation  Storage  (F.4S)  scheme  of  Brandt  [21] 
to  define  the  equivalent  fine  grid  problem  on  a  coarse  grid.  The  grid  transfer  operators 
for  the  solution,  residual,  and  coarse  grid  corrections  are  those  introduced  by  Jameson 
[1].  In  order  to  execute  the  multigrid  strategy  we  employ  a  fixed  W-type  cycle  with  two 
sweeps  on  each  coarse  grid.  To  provide  a  well  conditioned  starting  solution  for  the  fine 
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mesh  a  Full  Multigrid  (FMG)  method  is  used.  The  FMG  is  analogous  to  grid  sequencing; 
except  that  multigrid  cycles  are  performed  on  each  coarse  grid. 

Some  of  the  additional  elements  of  the  multigrid  method  are  not  necessarily  standard. 
A  smoothing  of  the  coarse  grid  corrections  being  transferred  to  the  finest  grid  was  found 
to  be  beneficial  in  transonic  computations  [17].  The  smoothing  was  accomplished  with 
the  implicit  residual  smoothing  mentioned  previously  and  a  constant  coefficient  ^  «  0.1. 
This  smoothing  of  the  residuals  on  the  way  to  finer  meshes  is  crucial  for  the  convergence 
of  the  multigrid  for  hypersonic  flows.  Such  a  process  acts  to  reduce  high  frequency  oscil¬ 
lations  caused  by  the  interpolation.  Hence,  it  becomes  especially  important  near  strong 
shocks,  where  nonphysical  upstream  influence  can  occur.  It  should  be  emphasized  that 
choosing  the  smoothing  parameter  too  large  can  slow  convergence.  Another  important 
element  for  hign  Mach  number  (j\/  >  10)  flows  is  the  coarse  grid  correction  scheme.  That 
is,  the  physical  viscous  terms  should  also  be  computed  on  the  coarse  meshes.  Difficulties 
with  applying  a  turbulence  model  on  a  coarse  grid  can  be  avoided  by  interpolating  the 
turbulent  viscosity  from  the  \'alues  on  the  finest  grid. 

Boundary  Conditions  and  Initialization 

-At  a  solid  surface  (wall)  boundar\'  the  no-slip  condition  is  enforced.  The  wall  pres¬ 
sure  is  set  to  the  v'alue  at  the  first  interior  solution  point,  and  thus,  a  reduced  normal 
momentum  equation  is  satisfied.  .An  adiabatic  wall  is  assumed.  In  a  finite-volume  for¬ 
mulation,  this  amounts  to  treating  the  Cartesian  velocity  components  as  antisymmetric 
functions  and  the  temperature  as  a  symmetric  function  with  respect  to  the  wall.  For  each 
of  the  physical  problems  considered  the  Mach  number  at  the  inflow  boundary  exceeds 
1.0.  Consequently,  the  dependent  variables  are  specified  at  this  boundary  according  to 
the  flow  conditions.  At  any  outflow  boundary,  we  apply  simple  extrapolation  of  the 
components  of  the  solution  vector.  In  general,  for  hypersonic  flows,  numerical  difficulties 
are  experienced  at  the  start  of  a  calculation  if  the  discrete  flow  field  is  initialized  with 
free-stream  conditions.  To  avoid  these  difficulties  we  apply  the  following  procedure.  The 
Mach  number  of  the  flow  is  set  to  a  lower  value  than  the  required  one.  Then  the  Mach 
number  is  gradually  increased  over  a  few  hundred  time  steps  until  the  desired  flow  con¬ 
ditions  are  obtained.  This  Mach  number  ramping  is  only  done  on  the  coarsest  mesh  in 
the  FMG  sequence. 

Results 

We  consider  the  following  two  test  cases:  (1)  two-dimensional  turbulent  flow  about 
the  N.AC.A  0012  airfoil,  and  (2)  three-dimensional  turbulent  flow  about  a  blunt  biconic. 
For  these  cases  we  found  a  significant  sensitivity  to  the  grid.  With  a  grid  suitable  for 
transonic  flows  (i.e.,  a  standard  C-type  mesh),  difficulties  in  convergence  were  e.xpe- 
riencod  at  high  Mach  numbers.  Furthermore,  the  usual  mesh  density  for  airfoil  flows 
yielded  a  poor  resolution  of  the  flow  field.  Tlie  grids  which  were  generated  for  these 
examples  all  follow  in  a  general  sense  the  bow  shock.  The  bow  shock  is  not  aligned  with 
the  grid,  but  as  a  minimum,  the  angle  between  the  bow  shock  and  the  grid  is  not  too 
large.  We  also  found  it  important  to  control  the  normal  stretching  rate  of  the  grid  in 
the  in  viscid  region  of  the  flow  field. 

For  the  first  case  we  consider  turbulent  flow  (a  Reynolds  number  Rcc  =  10")  about 
the  N.ACA  0012  airfoil  at  two  different  Mach  numbers.  In  all  the  calculations,  transition 
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was  assumed  to  occur  at  the  20  percent  chord  location.  This  ensured  that  the  algebraic 
turbulence  model  of  Baldwin  and  Lomax  {22j  did  not  determine  an  unrealistic  length 
scale  in  the  vicinity  of  the  leading  edge  due  to  the  presence  of  the  shock.  The  grids  for 
the  computations  were  generated  with  the  method  of  [23]  -  The  fine  grid,  consisting  of 
320  X  6-i  cells,  is  displayed  in  figure  1.  In  the  streamwise  direction,  the  grid  is  clustered  at 
the  leading  and  trailing  edges  of  the  airfoil,  with  a  tangential  spacing  of  approxi match' 
1  X  10~^  chords.  There  is  clustering  in  the  normal  direction  adjacent  to  the  surface,  so  as 
to  pro\  idc  adequate  resolution  of  the  turbulent  boundary  layer.  The  minimum  normal 
spacing  occurs  at  the  airfoil  leading  edge  and  is  1  x  lO"’  chords.  To  obtain  reasonable 
shock  definition  everywhere,  there  is  weak  normal  clustering  in  the  inviscid  region. 

Figures  2  and  3  show  the  Mach  number  and  pressure  contours,  respectively,  when 
—  10.  There  is  fairly  good  resolution  of  the  bow’  shock,  and  the  contours  are 
smooth.  The  pressure  contours  delineate  the  expected  fish  tail  shock  in  the  wake  region. 
Surface  pressure  and  skin- friction  distributions  for  this  case  are  presented  in  figure  4. 
In  addition,  velocity  profiles  at  the  midchord  and  9-5  percent  chord  locations  are  shown 
in  figure  -5.  The  same  results  computed  on  a  160  x  -32  mesh,  which  is  a  proper  subset 
of  the  320  X  64  mesh,  are  also  included  in  these  figures.  One  can  clearly  sec  that  there 
is  fairly  good  agreement  between  the  predictions  on  the  two  meshes.  Thus,  additional 
mesh  refinement  will  serve  priniariK  to  improve  the  downstream  resolution  of  the  bow- 
shock.  In  figure  6  the  \ariations  of  the  Mach  number  and  pressure  along  the  stagnation 
streamline  are  shown.  Comparing  the  pressure  jump  and  the  stagnation  pressure  with 
one-dimensional  theory,  we  find  the  differences  arc  less  than  2  percent.  With  the  cell- 
centered,  finite  volume  scheme,  one  should  keep  in  mind  that  the  computational  points 
do  not  lie  precisely  along  the  stagnation  streamline. 

Convergence  histories  for  the  Mach  10  case  arc  given  in  figure  7.  On  the  finest 
mesh  the  logarithm  of  the  residual  of  the  continuity  equation  is  reduced  7.5  orders  of 
magnitude  in  200  multigrid  cycles.  This  corresponds  to  an  average  rale  of  reduction  of 
0.917.  The  calculation  required  less  than  6  minutes  on  a  Cray  II  computer.  It  should 
be  emphasized  that  for  engineering  accuracy  (i.c.,  residual  reduced  by  3  orders)  the 
finest  mesh  calculation  required  less  than  3  minutes  of  CPU  time.  Note  that  when 
engineering  accuracy  is  achieved,  there  is  no  appreciable  improvement  in  the  viscous 
solution  accuracy  by  further  residual  reduction. 

In  figures  S-ll  results  for  =  20  and  Rcc  =  10~  are  presented.  Again  there  is  fairly 
good  representation  of  the  bow  shock,  and  the  Mach  and  pressure  contours  arc  smooth. 
Here  also  there  is  less  than  a  2  percent  difference  between  the  computed  pressure  jump 
and  the  stagnation  pressure  determined  along  the  stagnation  streamline  and  the  values 
from  one-dimensional  theory.  As  indicated  in  figure  12,  the  logarithm  of  the  residual 
is  reduced  nearly  5  orders  in  200  cj'clcs  on  the  finest  mesh.  There  is  a  slowdow-n  in 
the  asymptotic  convergence  rate  for  this  case.  At  this  point  the  precise  reason  for  such 
behavior  is  not  clear. 

We  now  consider  the  three-dimensional  case.  .A  mesh  consisting  of  128  x  96  x  24 
cells  was  used  for  the  blunt  biconic.  In  figure  13,  a  two-dimensional  slice  of  the  mesh  is 
depicted.  This  was  constructed  using  an  algebraic  mesh  generator.  A  tangent  hyperbolic 
distribution  function  was  used  in  the  normal  direction  to  obtain  an  acceptable  mesh  in 
both  the  viscous  and  inviscid  regions.  For  this  turbulent  flow  case,  =  6,  o  =  5 
degrees,  and  the  Revnolds  number  based  on  nose  diameter  is  2.S9  x  10*.  Figures  14  and 
15(a)  •  15(b)  show  Mach  number  and  pressure  contours,  respectively,  in  the  symmetry 


10 


plane.  One  sees  the  small  standoff  distance  of  the  shock  at  the  nose  and  how  the  shock 
remains  relativel)'  close  to  the  geometry  downstream.  The  shock  surface  is  delineated 
in  the  three-dimensional  view  of  pressure  contours  given  in  figure  15(c).  A  comparison 
is  made  in  figure  16  between  computed  pressure  distributions,  for  the  windward  and 
leeward  planes,  and  the  experimental  data  of  [24].  There  is  generally  good  agreement 
with  the  data.  In  figure  17  the  conx^ergence  history  is  presented  for  this  computation. 
The  logarithm  of  the  residual  has  been  reduced  about  5.5  orders  of  magnitude  in  about 
565  work  units,  which  roughly  corresponds  to  350  multigrid  cycles  on  the  finest  mesh. 
For  engineering  accuracy,  85  cycles  are  required  on  the  finest  mesh,  and  the  CPU  time 
is  about  30  minutes. 

Concluding  Remarks 

A  multigrid  method  with  central  differencing  has  been  successfully  applied  to  the 
solution  of  hypersonic  viscous  flows.  An  explicit  five  stage  R.unge-Kutta  scheme  has  been 
used  as  a  smoother  in  solving  the  time-dependent,  thin-layer  Navier-Stokes  equations. 
In  this  paper,  considerable  emphasis  has  been  focussed  on  the  dissipative  characteristics 
of  the  driving  scheme  for  the  multigrid  process.  The  presence  of  strong  shocks  has 
required  the  introduction  of  a  switching  function  for  the  numerical  dissipation  based 
on  TVD  principles.  This  nonlinear  switching  function  is  also  required  for  the  coarse 
grid  dissipation  model.  To^  importance  of  the  physical  diffusion  limit  on  the  time  step 
for  computing  hypersonic  flows  has  been  discussed.  By  introducing  appropriate  variable 
coefficients  for  implicit  residual  smoothing,  we  have  removed  the  diffusion  limit  and  have 
improved  the  convergence  significantly. 

Numerical  solutions  have  been  obtained  for  two-  and  three-dimensional  hypersonic 
turbulent  flows.  The  agreement  between  predictions  for  flow  over  a  blunt  biconic  compare 
well  with  experimental  data.  Engineering  accuracy  has  been  obtained  rapidly  in  all 
computations,  requiring  less  than  3  CPU  minutes  on  a  Cray  II  for  two-dimensional  cases 
and  about  30  CPU  minutes  for  a  three-dimensional  case. 
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(a)  Complete  airfoil 


(b)  Leading  edge  region 


Figure  Pressure  contours  for  Mach  10  turbulent 
flow  ovu  NACA  0012  airfoil  (a  =  0  deg,  Rcc  =  10^ 
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(b)  Surface  skin-friction  coefficients 


Figure  4  Surface  pressure  and  skin-friction  coefficients  for  Mach 
10  turbulent  flow  over  NACA  0012  airfoil  (a  =  0  cleg,  Rcc  =  10^) 
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(a)  X/C  =  0.5  (b)  X/C  =  0.95 

Figure  5  Velocity  profiles  for  Mach  10  turbulent 
flow  over  NACA  0012  airfoil  (a  =  0  deg,  Reg  =  10^) 


Figure  6  Pressure  and  Mach  number  variations 
along  stagnation  streamline  for  Mach  10  turbulent 
flow  (NACA  0012  airfoil,  a  =  0  deg,  Rcc  =  10^) 


Figure  7  Convergence  histories  for  Mach  10  turbulent 
flow  (NACA  0012  airfoil,  o  =  0  deg,  Rce  =  10^) 


(a)  Complete  airfoil 


(b)  Leading  edge  region 


Figure  8  Mach  number  contours  for  Mach  20  turbulent 
flow  over  NACA  0012  airfoil  (a  =  0  deg,  Rec  =  10^) 
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la)  Surface  pressure  coefficients  (b)  Surface  skin-friction  coefficients 

Figure  10  Surface  pressure  and  skin-friction  coefficients  for  Mach 
20  turbulent  flow  over  NACA  0012  airfoil  (a  =  0  deg,  Rcc  =  10^) 


Figure  11  Pressure  and  Mach  number  variations  Figure  12  Convergence  histories  for  Mach  20  turbulent 

along  stagnation  streamline  for  Mach  20  turbulent  flow  (NACA  0012  airfoil,  a  =  0  deg,  Rcc  =  10^) 

flow  (NACA  0012  airfoil,  a  =  0  deg,  TJce  =  10^) 
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(a)  Blowup  of  nose  region 


(b)  Symmetry  plane 


Figure  15  Pressure  contours  for  Mach  6  turbulent 
flow  over  blunt  biconic  (a  =  5  deg,  7iep=  2.89  x  10®) 
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(a)  Wndward  plane  (b)  Leeward  plane 

Figure  16  Pressure  coefficient  distributions  in  axial  direction  for  Mach 
6  turbulent  flow  over  blunt  biconic  (or  =  5  deg,  Rep=  2.89  x  10®) 
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Figure  17  Convergence  history  (blunt  biconic,  o  =  5  deg,  Rep  =-2.89  x  10®) 
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